function decisionValueDynamic()

[v, d] = decisionValueDynamic0(nb, r_p, r_n, r_s, nT,1);
figure;
plot((1:nT), d,'o-'); 
xlim([1 nT]);
ylim([0.5 1]);
xlabel('Time t');
ylabel('Policy \pi');
set(gcf,'paperunits','inches');
set(gcf,'papersize',[11 11]);
set(gcf,'paperposition',[0,0,11,11]);
end

function [optV optD] = decisionValueDynamic0(nb, r_p, r_n, r_s, nT, deltaT)
%************************************************************************
%POMDP VALUE FUNCTION FOR SEQUENTIAL DECISION PROBLEM
%To Infer The Hidden Binary State s = {+1, -1},
%From The Binary Observation State o = {+1, -1}
%With UNKOWN Emission Probability Model
%Assumption:
%     Belief b_{k+1}  = ( 1 - gamma_k) * b_k + gamma_k * 1_{o_k = s_k}
%            b_k is the belief at k step and
%            gamma_k is the learning rate at k step, inverse to time k. 
%Inputs:
%       r_p, reward for making the correct action
%       r_n, reward for making the incorrect action
%       r_s, reward for sampling
%Outputs:
%      optV, nb * nT time dependent optimal value function
%      optD, nT, time dependent optimal policy
%************************************************************************

if nargin < 6
    nb = 21;
    nT = 1;        %nT estimated points
    deltaT = 30;    %incremental intial time step
    r_p = 20;
    r_n = -100;
    r_s = -0.2;
end

b =  linspace(0,1,nb);
tol = 1.0e-4; %Transition matrix won't effectively change after maxT

%Reward function if no sampling is allowed, or the value function when the
%threshold is set to 0.5
R = max(b * r_p + (1 - b) * r_n, b * r_n + ( 1 - b) * r_p);

%Compute the value function for different policies
d = 1.0 * (nb - 2) / (nb - 1) / 2 + 1.0 * [1:(nb-1)/2] / (nb -1);

optD = zeros(1, nT);
optV = zeros(nb,nT);
%maxT = 10000;
% T = zeros(nb,nb, maxT);
% 
% for i = 1: maxT
%     T(:,:,i) =  transition(i, b);
% end

for initT = 1 : nT
    %For each starting time, find the best policy
    V = zeros(nb, (nb-1)/2);
    for d_i = 1 : (nb-1)/2 %for all possilbe policy
        m = true(nb,1);  %m(i) if i is a non-terminal state
        %state i is a terminal state if i is beyond the threshold set by the
        %policy
        m(b > d(d_i)) = false;
        m(b < 1 - d(d_i)) = false;
         %Receive reward if already in the terminal state
        T = transition(initT * deltaT, b);
        V(~m,d_i) =  R(~m);

        V(m, d_i) = (eye(sum(m)) - T(m,m)) \ (T(m,~m) * R(~m)' + r_s);
        
%         V(:,d_i) =  V(:,d_i) + r_s * (initT - 1);
%         t = (initT - 1 ) * deltaT; % Intitial time step
%         F = eye(nb);
%         while max(max(abs(F(m,m)))) > tol  %Value iteration              
%             t = t + 1;  
%             F(:,~m) = 0;
%             F(~m,:) = 0;
%             if t < maxT
%                 F(m,:) = F(m,:) * T(:,:,t);
%             else
%                 F(m,:) = F(m,:) * T(:,:,maxT);
%             end
%             V(m,d_i) = sum(F(m,m),2)* r_s +  V(m,d_i) + F(m,~m) * R(~m)'; 
%         end
    end

    [optV(:,initT), tempD] = max(V,[],2);
    optD(initT) = d(tempD((nb+1)/2));
    
%     if initT == 1 %print only when initT = 1
%         figure; 
%         plot(b,V);
%         hold on; 
%         plot(b,R, 'k', 'linewidth',4);
%     end
%         if initT == 10 %print only when initT = 1
%         figure; 
%         plot(b,V + 9 * r_s);
%         hold on; 
%         plot(b,R, 'k', 'linewidth',4);
%         end

    
end % End of initT

end %End of function

function T= transition(t,b)%Transition function should be time dependent
nb = length(b);
T = zeros(nb,nb);
for i = 1 : nb %P(s = S_R)
    %Probability that the next belief is b2R(next observation is O_R)?
    P_b2R = b(i);
    %Probability that the next belief is b2L(next observation is O_L)?
    P_b2L = 1 - P_b2R;    
    %NonBayesian belief update
    if P_b2R ~= 0
          b2R = beliefUpdate(b(i), t, 1); %The next possible belief when o_t+1 = O_R
    else
        b2R = 0;
    end
    if P_b2L ~=0
         b2L  = beliefUpdate(b(i), t, 0); %The next possible belief when o_t+1 = O_L
    else
        b2L = 0;
    end
    %Floors and ceils of b2R and b2L
    jR_F = floor(b2R * (nb-1)) + 1;
    jR_C = ceil(b2R * (nb-1)) + 1;
    jL_F = floor(b2L * (nb-1)) + 1;
    jL_C = ceil(b2L * (nb-1)) + 1;
    %Weights on each grid points
    a_R = b2R * (nb - 1) + 1 - jR_F;
    a_L = b2L * (nb - 1) + 1 - jL_F;
    T(i,jR_F) = T(i,jR_F) + P_b2R * (1 - a_R);
    T(i,jR_C) = T(i,jR_C) + P_b2R * a_R;
    
    T(i,jL_F) = T(i,jL_F) + P_b2L * (1 - a_L);
    T(i,jL_C) = T(i,jL_C) + P_b2L * a_L;
end
end

